---
title: "Main Observational Models"
---

# Setup

```{r packages and load data, message=F, warning=F}

knitr::opts_chunk$set(include = F)

## Packages
library(interplot)
library(stargazer)
library(sjPlot)
library(sjmisc)
library(tidyverse)
library(conflicted)
library(ggsci)
library(ggpubr)
library(knitr)
library(broom.mixed)
library(ggeffects)
library(marginaleffects)
library(haven)
library(here)
library(tictoc)
library(glmmTMB)
library(texreg)
library(dotwhisker)
library(parallel)
library(survey)
library(interplot)

# Make things a bit faster
cores <- detectCores() - 2


conflict_prefer("select", "dplyr")
conflict_prefer("filter", "dplyr")

# Check directory
here()

```


# Load Data and Minor Recodes

```{r load data and make factors}


## Load Data 

load(here("data/lapop_10_18.Rdata")) # slim data -- recodes only


# DATA PREP 

# Select Variables
vars_close_cong <- c(
  "close.cong", "dem.best", "sat_dem", "winner", "non_voter", "pres_support",
  "edu", "wealth", "rural", "wealth", "female", "trust.people", "neigh.safe",
  "econ.eval", "dem.best.tj", "pop.speech", "maj", "v2x_libdem",
  "mass.polr", "country", "year", "year", "strata", "weight1500"
)


vars_close_court <- c(
  "close.court", "dem.best", "sat_dem", "winner", "non_voter", "pres_support",
  "edu", "wealth", "rural", "wealth", "female", "trust.people", "neigh.safe",
  "econ.eval", "dem.best.tj", "pop.speech", "maj", "v2x_libdem",
  "mass.polr", "country", "year", "year", "strata", "weight1500"
)


dichot.vars.cong <- c("close.cong", "rural", "female", "non_voter", "winner", "year")
dichot.vars.court <- c("close.court", "rural", "female", "non_voter", "winner", "year")


# Select vars and make factors
cong_dat <- lapop.cleaned %>%
  select(any_of(vars_close_cong)) %>%
  na.omit() %>% 
  mutate(across(all_of(dichot.vars.cong), factor)) %>%
  zap_formats() %>% 
  zap_labels()

court_dat <- lapop.cleaned %>%
  select(any_of(vars_close_court)) %>%
  na.omit() %>% 
  mutate(across(all_of(dichot.vars.court), factor)) %>% 
  zap_formats() %>% 
  zap_labels()



```

# Close Congress Models

```{r close congress models}
## No interaction 

cong_mod1 <- as.formula("close.cong ~ dem.best + sat_dem + edu + wealth + rural + wealth + female + 
                      trust.people + neigh.safe + econ.eval + winner + non_voter + pop.speech +
                      maj + v2x_libdem + mass.polr + dem.best.tj + (1 | country:year) + 
                      (winner + non_voter | country)")


cong_fit1 <- glmmTMB(formula = cong_mod1,
                     family = binomial(link = "logit"), 
                     data = cong_dat, 
                     control=glmmTMBControl(parallel = cores))


## Supporter x Populist Speech

cong_mod2 <- as.formula("close.cong ~ dem.best + sat_dem + edu + wealth + rural + wealth + female + 
                      trust.people + neigh.safe + econ.eval + non_voter + pop.speech * winner +
                      maj + v2x_libdem +  mass.polr + dem.best.tj + (1 | country:year) +
                      (winner + non_voter | country) + year")


cong_fit2 <- glmmTMB(formula = cong_mod2,
                     family = binomial(link = "logit"), 
                     data = cong_dat, 
                     control=glmmTMBControl(parallel = cores))



## Supporter x Pop Speech + Non-Voter x Pop Speech 

cong_mod3 <- as.formula("close.cong ~ dem.best + sat_dem + edu + wealth + rural + wealth + female + 
                      trust.people + neigh.safe + econ.eval +  pop.speech * winner + pop.speech * non_voter + 
                      maj + v2x_libdem + mass.polr + dem.best.tj + (1 | country:year) +
                      (winner + non_voter | country) + year")


cong_fit3 <- glmmTMB(formula = cong_mod3,
                     family = binomial(link = "logit"), 
                     data = cong_dat, 
                     control=glmmTMBControl(parallel = cores))


cong_mods1 <- list(cong_fit1, cong_fit2, cong_fit3)
```

## Regression Tables

```{r close congress models}

screenreg(cong_mods1, 
             stars = c(.1, .05, .01))

coef_labs_main_si <- list("pop.speech" = "Populist Discourse", "pop.speech:winner1" = "Populist Speech $\\times$ Pres. Supporter",
             "pop.speech:non_voter1" = "Populist Speech $\\times$ Non-Voter", "winner1" = "Presidential Supporter", 
             "non_voter1" = "Non_Voter","dem.best" = "Support for Democracy", "sat_dem" = "Satisfaction with Democracy", 
             "edu" = "Education", "wealth" = "Wealth", "rural1" = "Rural", "female1" = "Female", 
             "trust.people" = "Interpersonal Trust", "neigh.safe" = "Neighborhood Safe", 
             "econ.eval" = "Economic Evaluation", "maj" = "Legislative Control", 
             "v2x_libdem" = "Liberal Democracy Index", "mass.polr" = "Mass Polarization", 
             "dem.best.tj" = "Mass Support for Democracy")

coef_labs_main <- list("pop.speech" = "Populist Discourse", "pop.speech:winner1" = "Populist Speech $\\times$ Pres. Supporter",
             "pop.speech:non_voter1" = "Populist Speech $\\times$ Non-Voter", "winner1" = "Presidential Supporter",
             "non_voter1" = "Non_Voter")



screenreg(cong_mods1, file = here("tables/main_cong_mods.txt"))

texreg(cong_mods1,
       custom.coef.map = coef_labs_main_si,
       caption = "Multilevel logistic models of support for closing congress. Coefficients are in logged odds.",
       caption.above = T,
       label = "tab:cong_mods1_si",
       stars = c(.1, .05, .01),
       file = here("tables/main_cong_mods_si.tex"))

texreg(cong_mods1,
       custom.coef.map = coef_labs_main,
       caption = "Multilevel logistic models of support for closing congress. Coefficients are in logged odds.",
       caption.above = T,
       label = "tab:cong_mods1",
       stars = c(.1, .05, .01),
       file = here("tables/main_cong_mods.tex"))


```


# Close Court Models

```{r close court models}
## No interaction 
court_mod1 <- as.formula("close.court ~ dem.best + sat_dem + edu + wealth + rural + wealth + female +
                        trust.people + neigh.safe + econ.eval + pop.speech + winner + non_voter +
                        maj + v2x_libdem +  mass.polr + dem.best.tj + (1 | country:year) +  
                        (winner + non_voter | country) + year")

court_fit1 <- glmmTMB(
  formula = court_mod1,
  data = court_dat,
  family = binomial(link = "logit"))


## Populism x Presidential Supporter 
court_mod2 <- as.formula("close.court ~ dem.best + sat_dem + edu + wealth + rural + wealth + female +
                        trust.people + neigh.safe + econ.eval + non_voter + pop.speech * winner +
                         maj + v2x_libdem +  mass.polr + dem.best.tj + (1 | country:year) + 
                        (winner + non_voter | country) + year")

court_fit2 <- glmmTMB(
  formula = court_mod2,
  data = court_dat,
  control=glmmTMBControl(parallel = cores),
  family = binomial(link = "logit"))


##  Supporter x Pop Speech + Non-Voter x Pop Speech  
court_mod3 <- as.formula("close.court ~ dem.best + sat_dem + edu + wealth + rural + wealth + female +
                        trust.people + neigh.safe + econ.eval + pop.speech * winner + pop.speech * non_voter +
                         maj + v2x_libdem + mass.polr + dem.best.tj + (1 | country:year) +
                        (winner + non_voter | country) + year")

court_fit3 <- glmmTMB(
  formula = court_mod3,
  data = court_dat,
  control=glmmTMBControl(parallel = cores),
  family = binomial(link = "logit"))



## Models Table S4
court_mods1 <- list(court_fit1, court_fit2, court_fit3)
```

## Regression Tables

```{r close court models}
# Regression Tables

screenreg(court_mods1)

screenreg(court_mods1, file = here("tables/main_court_mods.txt"))

texreg(court_mods1,
       custom.coef.map = coef_labs_main_si,
       caption = "Multilevel logistic models of support for the court. Coefficients are in logged odds.",
       caption.above = T,
       label = "tab:courtmods1_si",
       stars = c(.1, .05, .01),
       file = here("tables/main_court_mods_si.tex"))

texreg(court_mods1,
       custom.coef.map = coef_labs_main,
       caption = "Multilevel logistic models of support for closing the court. Coefficients are in logged odds.",
       caption.above = T,
       label = "tab:courtmods1",
       stars = c(.1, .05, .01),
       file = here("tables/main_court_mods.tex"))

```


# Plot Predictions

## Plot Close Cong

```{r close cong model predictions}


## Get predictions 
cong_preds <-  predictions(cong_fit3, 
            newdata = datagrid(winner = c("0", "1"),
                                non_voter = c("0", "1"),
                                pop.speech = seq(0, 1.75, .05)))


cong_preds <- cong_preds %>% 
  select(estimate, conf.low, conf.high, winner, non_voter, pop.speech)

# Label for plotting 
cong_preds <- cong_preds %>% 
  mutate(group = case_when(winner == 0 & non_voter == 1 ~ "Non-Voter",
                           winner == 1 & non_voter == 0 ~ "Presidential Supporter",
                           winner == 0 & non_voter == 0 ~ "Opposition Supporter"))

cong_preds$group <- factor(cong_preds$group, levels = c("Opposition Supporter",
                                                        "Non-Voter",
                                                        "Presidential Supporter"))

# drop winner = 1, non-voter = 1
cong_preds <- na.omit(cong_preds)

# plot:
cong_plot <- ggplot(cong_preds, aes(x = pop.speech, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_line() +
  geom_ribbon(alpha = .2) +
  scale_x_continuous(limits = c(0, 1.75), breaks = seq(0,1.75, .5)) +
  facet_wrap(~group) + 
  labs(title = "Justified to Close Congress") +
  theme_pubclean() +
  theme(plot.title = element_text(hjust = 0.5), 
        axis.title = element_blank())


```

## Plot Close Court

```{r close court model predictions}


# get predictions
court_preds <-  predictions(court_fit3, 
                       newdata = datagrid(winner = c("0", "1"),
                                          non_voter = c("0", "1"),
                                          pop.speech = seq(0, 1.75, .05)))


court_preds <- court_preds %>% 
  select(estimate, conf.low, conf.high, winner, non_voter, pop.speech)

# label for plotting
court_preds <- court_preds %>% 
  mutate(group = case_when(winner == 0 & non_voter == 1 ~ "Non-Voter",
                           winner == 1 & non_voter == 0 ~ "Presidential Supporter",
                           winner == 0 & non_voter == 0 ~ "Opposition Supporter"))


court_preds <-  na.omit(court_preds)

court_preds$group <- factor(court_preds$group, levels = c("Opposition Supporter",
                                                        "Non-Voter",
                                                        "Presidential Supporter"))


court_plot <- ggplot(court_preds, aes(x = pop.speech, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_line() +
  geom_ribbon(alpha = .2) +
  facet_wrap(~group) + 
  labs(title = "Justified to Close Court") +
  scale_x_continuous(limits = c(0, 1.75), breaks = seq(0,1.75, .5)) +
  theme_pubclean() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.title = element_blank())


pred_probs <- ggarrange(cong_plot, court_plot, nrow = 2)
```


## Combine and Save Figs

```{r save predictions}

# annotate
annotate_figure(pred_probs, 
                bottom = text_grob("Populist Discourse"), 
                left = text_grob("Probality Justified", rot = 90)
                )


# Save to dropbox
ggsave(
  here("figures/all_preds.pdf"), 
  height = 6, 
  width = 8)


```



# Discussed in Text

## Summary Stats

```{r  discussed in text}
# Summary Stats

# Means of close court and close cong by country year
cong_means <- cong_dat %>%   
  mutate(close_cong_num =  
           case_match(close.cong, "1" ~ 1, "0" ~ 0)) %>% 
  group_by(country, year) %>% 
  summarise(mean_close_cong = mean(na.omit(close_cong_num))) %>% 
  arrange(mean_close_cong)

# Means by country year
court_means <- court_dat %>% 
  mutate(close.court.num =  
           case_match(close.court, "1" ~ 1, "0" ~ 0)) %>% 
  group_by(country, year) %>% 
  summarise(mean_close_court = mean(na.omit(close.court.num))) %>% 
  arrange(mean_close_court) %>% 
  arrange(mean_close_court)


# Mean support closing congress
cong_dat %>% 
  select(close.cong) %>% 
  na.omit() %>% 
  group_by(close.cong) %>% 
  summarise(n = n()) %>% 
  mutate(prop = n/sum(n))

# Mean support closing court
court_dat %>% 
  select(close.court) %>% 
  na.omit() %>% 
  group_by(close.court) %>% 
  summarise(n = n()) %>% 
  mutate(prop = n/sum(n))

# Summary stats for populist discourse
cong_dat %>% 
  select(country, year, pop.speech) %>% 
  unique.data.frame() %>% 
  summarise(mean = mean(pop.speech),
            sd = sd(pop.speech), 
            min = min(pop.speech), 
            max = max(pop.speech))


```


## Predictions

```{r  discussed in text}
## Effect of moving from min to max using marginaleffects package 

f_diff_cong <-  avg_comparisons(cong_fit1, variables = list(pop.speech = c(0, 1.75))) 

f_diff_court <- avg_comparisons(court_fit1, variables = list(pop.speech = c(0, 1.75))) 

## Check some summary stats 
cong_dat_country_yr <- cong_dat %>% 
  select(pop.speech, country, year) %>% 
  unique.data.frame() 

court_dat_country_yr <- court_dat %>% 
  select(pop.speech, country, year) %>% 
  na.omit() %>% 
  unique.data.frame() 

summary(cong_dat_country_yr$pop.speech)
summary(court_dat_country_yr$pop.speech)

# Predictions Close Cong
discuss_cong <- predictions(cong_fit3, 
            newdata = datagrid(winner = c("0", "1"),
                                non_voter = c("0", "1"),
                                pop.speech = c(0, 1.75))) %>% 
  select(estimate, winner, non_voter, pop.speech) %>% 
  mutate(group = case_when(winner == 0 & non_voter == 1 ~ "Non-Voter",
                           winner == 1 & non_voter == 0 ~ "Presidential Supporter",
                           winner == 0 & non_voter == 0 ~ "Opposition Supporter")) 

discuss_cong %>% 
  arrange(pop.speech) %>% 
  na.omit() %>% 
  kable()

# Predictions Close Court
discuss_court <- predictions(court_fit3, 
            newdata = datagrid(winner = c("0", "1"),
                                non_voter = c("0", "1"),
                                pop.speech = c(0, 1.75))) %>% 
  select(estimate, winner, non_voter, pop.speech) %>% 
  mutate(group = case_when(winner == 0 & non_voter == 1 ~ "Non-Voter",
                           winner == 1 & non_voter == 0 ~ "Presidential Supporter",
                           winner == 0 & non_voter == 0 ~ "Opposition Supporter")) 

discuss_court %>% 
  arrange(pop.speech) %>% 
  na.omit() %>% 
  kable()





```




